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By means of molecular dynamics, we study the structure and the dynamics of a microscopic 
model for colloidal gels at low volume fractions. The presence of directional interactions leads to 
the formation of a persistent interconnected network at temperatures where phase separation does 
not occur. We find that large scale spatial correlations strongly depend on the volume fraction and 
characterize the formation of the persistent network. We observe a pre-peak in the static structure 
factor and relate it to the network structure. The slow dynamics at gelation is characterized by the 
coexistence of fast collective motion of the mobile parts of the network structure (chains) with large 
scale rearrangements producing stretched exponential relaxations. We show that, once the network 
is sufficiently persistent, it induces slow, cooperative processes related to the network nodes. We 
suggest that this peculiar glassy dynamics is a hallmark of the physics of colloidal gels at low volume 
fractions. 

PACS numbers: 82.70.Gg, 82.70.Dd, 64.70.Pf 



I. INTRODUCTION 

Soft solids obtained from suspensions of attractive col- 
loidal particles are of central relevance in fields ranging 
from food processing to the design of new smart materi- 
als . In most cases of practical interest these solids are 
not crystalline but gels or soft glasses [2] : they arise from 
the structural arrest of the suspension into a disordered 
state, in which the particles can support external stress 
and therefore behave macroscopically like a solid. The 
mechanical and rheological properties of this soft mat- 
ter depend crucially on its structural features which, in 
turn, are the result of a delicate interplay between the 
underlying thermodynamics and the arrest conditions. 
Hence, the understanding of the mechanisms leading to 
structural arrest is decisive for many technological ap- 
plications in which colloidal suspensions are involved. 
Whereas it is clear that in dense suspensions the crowd- 
ing of the particles is the main mechanism leading to 
structural arrest and yielding, suspensions of attractive 
colloids allow to form soft solids also at low volume frac- 
tions, where aggregation produces large scale structures 
P-fn]. At intermediate volume fractions (say above 10- 
15% depending on the specific system) possible scenarios 
for the arrest are the crowding of the aggregates act- 
ing as the new units of a dense glassy system [TJHIS], 
or the geometrically frustrated arrangements of locally 
stable large scale structures (such as lamellae or meso- 
scopic elongated aggregates) p!6UT9] . At sufficiently low 
volume fractions, instead, the structure of the arrested 
states, although created via reversible aggregation, re- 
sembles more the open fractal networks typically created 
by diffusion-limited aggregation processes [20 . In this 
situation, where the arrest seems to be intimately related 
to the network structure, the understanding of the arrest 
phenomenon is far from being reached. 



Our work addresses this last situation, which has re- 
ceived only a limited attention in past theoretical and 
numerical studies. Although experiments show that re- 
versible aggregation processes lead to the formation of 
open, persistent network structures, most of the estab- 
lished models for the attractive effective interactions, 
very successful for many aspects [21], fail to produce at 
these low volume fractions network structures that are 
stable against macro- or micro-phase separation. Con- 
focal microscopy images obtained in recent experiments 
[22H24] demonstrate that the gel networks formed at 
low volume fractions have a distribution of the particle 
coordination number n that is strongly peaked around 
n c^2^3, suggesting that the effective interactions include 
also many-body terms. 

As a matter of fact, in these systems there are several 
possible sources for such anisotropics, since the particle 
surface may not be smooth or the building blocks of the 
gel are not the primary particles but larger aggregates 
of irregular shape [11]. Because, up to now, there is no 
possible quantitative, microscopic derivation of effective 
interaction which accounts for these complications, we 
have recently proposed a phenomenological model con- 
taining the relevant qualitative ingredients. Following 
up the idea that directional terms in the effective inter- 
actions can lead to the formation of an open network 
that is thermodynamically stable, our phenomenological 
model for the effective interactions contains a directional 
term which introduces a local rigidity [25H27) . Our aim is 
not to reproduce one specific experimental system but to 
understand how this type of effective interaction might 
be relevant for the gelation phenomena as observed in 
diluted suspensions of attractive colloidal particles. 

In our model, without imposing a maximal local con- 
nectivity [28] -[30], the network structure arises from the 
balance between the internal energy and the entropy, de- 
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pending on temperature and volume fraction. In a first 
study we have shown [26l [27] that the formation of the 
persistent network produces the coexistence, in the gel, 
of very different relaxation processes at different length 
scales: the relaxation at high wave vectors is due to the 
fast cooperative motion of pieces of the gel structure (e.g. 
the chains connecting two nodes), whereas at low wave 
vectors the overall rearrangements of the heterogeneous 
gel make the system relax via a stretched exponential 
decay of the time correlators. The coexistence of such 
diverse relaxation mechanisms is characterized by a typ- 
ical crossover length which is of the order of the network 
mesh size [31j. 

In this paper we present a systematic molecular dy- 
namics study of the model for different temperatures and 
volume fractions. We will argue that the structural and 
dynamical features of this model should be quite gen- 
eral and could be relevant to different experimental gel- 
forming systems. 

In section |ll| we give the details of the model and of 
the numerical simulations. The structure of the system 



is analyzed in section III and in section IV we study the 
dynamics. All the results are gathered into a coherent 
picture in section |V| 



II. MODEL AND NUMERICAL SIMULATIONS 

We consider identical particles of radius a, interacting 
via a potential K//- As mentioned above we seek a phe- 
nomenological model for the effective interactions which 
can lead to gelation in attractive colloidal suspensions at 
low volume fractions. For this, for the inter-particle at- 
traction we consider a Lennard- Jones-like potential of the 
form VLjiUj) = 2?>e{{(jLj/rijY^ - {(Tlj /rijY^), produc- 
ing a narrow attractive well (r^j is the distance between 
the centers of particles i and j). To introduce direction- 
ality, we have considered an additional soft-sphere repul- 
sion modulated by a geometric term, which reduces the 
attractive well of VLj(^ij) unless particles i and j are 
in specific relative configurations. To define such con- 
figurations the particles are decorated with sticky points 
occupying the vertices of the icosahedron inscribed in the 
sphere of diameter <j and centered in the particle posi- 
tion (see Figjl]). Between particles i and j the directional 
effect is then obtained in the following way: 
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(1) 



Here is the position of i, r^p gives the position of the 
point p on particle i. The sum runs over the 12 points of 

respectively i and j. f{ri — rjp) = that, if 

there is one of these points whose distance from the cen- 
ter of particle i is small compared to Vd{rij) = and 




FIG. 1: The cartoon illustrates the mechanism of the direc- 
tional interactions: 1) Two particles close within the attrac- 
tion range and touching with the sticky points correspond to 
an internal energy per particle of —1 (in units of ksT). Any 
other configuration corresponds to an internal energy per par- 
ticle of —0.2. 2) The angle a among 3 neighboring particle is 
unlikely to be smaller than 70 degrees. 



the additional soft sphere repulsion does not contribute. 
If this distance is larger, then Vd{rij) = f (^)"^^, reduc- 
ing the depth of the original attractive well to 20%. As 
a consequence, two bonded particles will be subjected to 
the full attraction Vlj only if they touch in those spe- 
cific configurations and a reduced attraction otherwise. 
The choice of 12 points allows us to define such relative 
configurations without imposing a local maximum con- 
nectivity to the particles: for example, the choice of 3 or 
4 points would anyway lead to the formation of an open 
network interfering with the effect of directionality which 
is more appropriate for the physics of these systems. 

The term Vd alone is not able to effectively limit the 
functionality of the particles at the volume fraction and 
temperatures considered here [33] and will lead to a gel 
network whose structure and dynamics are strongly af- 
fected by phase separation kinetics. Therefore, in our 
effective interactions we consider a third term which im- 
poses a certain angular rigidity to the bonds ij and ik: 



V2^{rijk) = 13. 5e 



2 \ 9 _ r^ (rfc-r,)-(r^-r^) 



(2) 

making that the angle among the three neighbor particles 
is unlikely to be smaller than a certain value a (see Figjl]). 
The resulting interaction potential is 



Veff = VLJ + Vd + Vs. 



(3) 



The choice of including both the terms Vd and Vs has 
been made in the spirit of investigating more deeply their 
relative contribution to the formation of the open net- 
work (and to its dynamics), at a stage where the presence 
itself of any directional effect was questioned in the sci- 
entific community and therefore not considered in most 
of the cases. Our investigation of the different role of 
the two terms, albeit partial, indicates that a simplified 
version of the model with only Vlj + Vs should lead to 
very similar structural and dynamical features [34] , 

We have implemented Vef / in a constrained molecular 
dynamics code to perform micro-canonical simulations. 
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using a suitable combination of the algorithms RATTLE 
and SHAKE ^35. 

Although the model may seem quite complicate, in the 
following sections it will be shown how it undoubtedly 
picks relevant aspects of the physics it was meant to de- 
scribe and it does contribute to its understanding. 

In Eqsjl][3) the parameters a^j, a, 6, d and a can be 
used to tune the effective interactions thus giving rise to a 
large variety of interesting possibilities (e.g. varying the 
connectivity from chains to compact structures, changing 
the local rigidity of the network...). However the desired 
features of the persistent structure do not correspond to 
a very specific combination of the parameters, but to 
reasonable ranges of values which allow for extensive in- 
vestigations. In particular, we have selected the ranges 
which allow to give rise to the desired features, like a 
relatively narrow attractive well, a connectivity greater 
than 2 and a certain angular rigidity. Here we consider 
one reasonable choice which guarantees the formation of 
a persistent open network, i.e. ctlj = 0.922, a = 1.0, 
d = 0.43, b = 0.34 and a = OAtt. We have already 
shown that with these parameters at a volume fraction 
(j) = 0.05, the static structure of this system at low tem- 
perature is an open network structure [25 , in qualitative 
agreement with the one of colloidal gels. Furthermore 
the system does not show any sign of phase separation 
in the temperature range investigated and its relaxation 
time increases rapidly with decreasing T, i.e. the static 
and dynamic properties of the system are indeed very 
similar to the ones found in real colloidal gels [6 . 

In the following sections we will discuss the results of 
molecular dynamics simulations in the micro-canonical 
ensemble, using a time step of 0.002. The unit of time 
is Y^mcr^/e, with m the mass of a particle. All the data 
refers to simulations performed with 8000 particles in cu- 
bic boxes of size L = 37.64, 43.09, 55.10 in unit of a, cor- 
responding respectively to particle densities of p = 0.15, 
0.1, and 0.05 from which we have estimated approxi- 
mately a volume fraction (j) 0.075, 0.05, and 0.025. For 
every value of L we have studied the system at the tem- 
peratures 1.0, 0.7, 0.5, 0.3, 0.2, 0.15, 0.1, 0.09, 0.08, 0.06, 

0. 055, and 0.05. In order to do this we have performed 
the following equilibration protocol: starting from ini- 
tial high temperature random configurations, the system 
is equilibrated at each temperature by replacing all the 
velocities of the center of mass of the particles and the 
angular velocities with values extracted from a Maxwell- 
Boltzmann distribution every A time steps (A is suit- 
ably varied with temperature from 10 to 10^ MD steps). 
We have checked that after equilibration the energy is 
constant, showing no significant drift over the simulation 
time window. We have monitored independently transla- 
tional and rotational contributions to kinetic energy. Fi- 
nally, we have verified that different one- and two- time 
autocorrelation functions reach the equilibrium behavior, 

1. e. do not age. 

From these equilibrated configurations we start the 
data production. The equilibration time grows accord- 



ingly to the relaxation time in the system. At the lowest 
temperatures the equilibration procedure required up to 
2-10^ MD steps. For each temperature and volume frac- 
tion we generated five independent runs, over which the 
results presented here have been averaged. 

III. STRUCTURE 

We investigate the structural changes in the system 
at the different temperatures by means of a comparative 
study of the density fiuctuations, of the local connectivity 
and of the aggregation process. 

A. Static structure factor 

After equilibrating the system at different tempera- 
tures, we calculate the static structure factor: 

jk 

where the values of the modulus of the wave vector q con- 
sidered are the ones compatible with the periodic bound- 
ary conditions of the simulations box. 

In Fig. 2, we plot S{q) as a function of q for differ- 
ent temperatures, at the volume fractions (j) = 0.025 and 
(j) = 0.075. The plots show that, upon lowering the tem- 
perature, major structural changes occur due to aggre- 
gation, as indicated by the peak located at a wave vector 
approximately corresponding to the first neighbor dis- 
tance within the potential well i.e. q 7.6, and, more 
importantly, by the simultaneous onset of spatial cor- 
relations at low wave vectors. At high wave vectors the 
change of S{q) from high to low temperatures looks qual- 
itatively the same at different 0. The spatial correlations 
at low wave vectors, and therefore the mesoscopic and 
large scale structures, seem instead to depend strongly 
not only on T but also on (j). We observe that the data of 
Fig. 2 show no sign of a phase separation, which would be 
indicated by a rapidly growing peak at low wave vectors. 
Instead we see at low q the onset of spatial correlations 
extending over different length scales, suggesting the for- 
mation of an extended disordered structure. We compare 
in Fig. [sjthe S{q) at the different (j): at high temperature 
(inset of the figure), spatial correlations in the system are 
mainly temperature-controlled, i.e. depend only weakly 
on (j). At the lowest temperature (main frame), instead, 
space correlations are dominated by temperature up to 
length scales of the order of a few particle diameters. In 
this region of i.e. 2.0 < < 5.0, S{q) shows the same 
(X 1/q dependence for the three values of volume frac- 
tions considered. This clearly indicates the presence of 
mesoscopic elongated structures (chains) as a result of 
the aggregation driven by the potential energy. Finally, 
beyond length scales of the order of 3-4 particle diameters 
(i.e. q < 2.0), the strength and range of spatial correla- 
tions are strongly controlled by the volume fraction. 
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FIG. 2: The static structure factor as a function of the wave 
vector q at volume fraction (j) = 0.025 (a) and (j) = 0.075 
(b). In each plot, from bottom to top (at low q), T = 
1.0, 0.5, 0.2, 0.1, 0.09, 0.08, 0.07, 0.06, 0.055, 0.05. The dashed 
line indicates the dependence 1/q. 



If we consider now S{q) in analogy with the static 
structure factor of a polymer chain solution [36], we can 
interpret the mesoscopic length scales 2.0 < g < 7.0 as an 
intra-molecular regime for spatial correlations due to the 
elongated aggregates and the macroscopic scales q < 2.0 
(i.e. for lengths up to the simulations box linear size) 
as an inter-molecular regime. Following this description, 
the intra-molecular regime is strongly controlled by the 
interaction potential (and therefore by temperature), and 
extends up to a length-scale which gives a rough estimate 
of the persistence length of the elongated aggregates (> 
3-4 particle diameters). Beyond the persistence length, 
the S{q) captures the inter-molecular regime: the fact 
that this regime depends quite dramatically on the vol- 
ume fraction even in the small interval here considered 
(A0 = 0.05) suggests that the linear size of the aggre- 
gates extends well beyond their persistence length. 



To conclude this section, we make two remarks. First, 
we notice that at the lowest volume fraction (j) = 0.025, 
one recognizes a power law regime in S{q) at intermedi- 
ate and low q. It is clear that this cannot be interpreted 
as a fractal regime, since the exponent would be less than 
1.0. It should instead correspond to the crossover from 
the inter-molecular regime to an eventual fractal regime 
which is only detectable for a larger system. This will 
be confirmed by the analysis of the structure in the fol- 
lowing sections. Second, we observe that, at the highest 
volume fraction, (j) = 0.075, S{q) displays at low T a well 
defined pre-peak at a wave vector q 2.0. In network- 
forming ionic liquids, this kind of pattern in scattering 
intensity is typically associated to the presence of a net- 
work structure [37 . In gelling colloidal suspensions it 
has also been related to the presence of stable clusters 
of a typical size of 3-6 particles [6l |38] and its connec- 
tion to gelation is debated. In fact, S{q) can give in- 
formation only on the typical distances between pairs of 
particles occurring in the structure and does not allow 
on its own to discriminate between pairs of particles in 
the same aggregate or from different aggregates. Here, 
the pre-peak arises from the crossover between the intra- 
molecular and the inter-molecular regime at this volume 
fraction due to two competing effects: pairs of particles 
in the same large aggregate (chain) will be typically sep- 
arated at least by a distance of 3-4 particle diameters, 
due to the local rigidity; pairs of particles separated by 
larger distances, belonging or not to the same aggregate, 
will be limited instead upon increasing the volume frac- 
tion. With numerical simulations we have the possibility 
to access further structural features, to complement the 
information obtained through S{q). Therefore in the fol- 
lowing we exploit it by analyzing the local connectivity 
of particles, as well as the cluster size distribution, as a 
function of the temperature and of the volume fraction. 



B. Connectivity 

We define the coordination number c(n) as the fraction 
of particles that have exactly n neighbors. Two particles 
are considered to be neighbors if their distance is less that 
^min=l-l5 the location of the first minimum in the radial 
distribution function. By means of c(n) we characterize 
the change in the topology of the structure occurring if 
T is decreased. In Fig. [4j c(n) is plotted as a function of 
inverse temperature for different values of n and for the 
different volume fractions considered here. 

As a general result, at high temperatures the vast ma- 
jority of the particles are isolated, (n = not shown), and 
the fraction of dimers, n = 1, is relatively high (20—30%). 
With decreasing T this fraction initially increases, but 
for T < 0.2 it starts to rapidly decrease with decreasing 
T. Our data shows that this happens because particles 
start to form larger structure which are mainly chains: 
the fraction of particles with exactly two nearest neigh- 
bors increases rapidly and these (local) configurations be- 
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FIG. 3: 5'(^) as a function of the wave vector at volume 
fractions = 0.025, 0.05, and 0.075 at T = 0.05 (main frame) 
and T — 1.0 (inset). The full line indicate the dependence 
l/q. 



come by far the most prevalent ones at low T. Last but 
not least also the fraction of particles with n = 3 neigh- 
bors increases with decreasing T. Therefore the curves of 




1/T 

FIG. 4: Coordination number c{n) as a function of inverse 
temperature. The different symbols correspond to n = 1 (o), 
n — 2 (□) and n = 3 (A). The different connecting lines 
correspond to = 0.025 (dashed), (j) = 0.05 (dotted) and 
= 0.075 (full). 




FIG. 5: Cluster size distribution n{s) for the different tem- 
peratures at volume fraction = 0.025 (a), and (j) = 0.075 
(b). The dashed line corresponds to s~^"^. 



cally increases with (j): this suggests that the formation 
of chains is strongly driven from the interaction, whereas 
the increasing steric hindrances between large aggregates 
due to increasing (j) will favor the formation of junctions 
between the chains by counterbalancing the angular re- 
pulsion Vz- 

C. Clusters 



Fig. |4] indicate that upon decreasing T the particles form 
chains which are connected by bridging points or nodes 
(n = 3) to form an open network. It is interesting to re- 
mark that at the lowest temperature there are practically 
no free particles. The fraction of dimers (n = 1) is minor 
and it significantly decreases with increasing the volume 
fraction. The fraction of particles with n = 2 does not 
strongly depend on (j) in the range of values explored, 
whereas the fraction of particles with n = 3 monotoni- 



We have monitored the aggregation process and char- 
acterized the structure of the system on large length 
scales using n{s)^ the number of clusters that have exactly 
s particles. We define that a particle belongs to a cluster 
if its distance from at least one member of the cluster is 
less than rmin, the location of the first minimum in the 
radial distribution function. This distribution is shown 
in Fig. 6 for all temperatures and volume fractions inves- 
tigated. For the different volume fractions, at high tem- 
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peratures (T > 0.3) the distribution nicely follows an ex- 
ponential law, a behavior that corresponds to the random 
formation of transient clusters of non-bonded particles at 
low densities. Around T = 0.1 the shape of the distribu- 
tion starts to change strongly, indicating that dimers are 
the most probable clusters. At the same time, n{s) also 
starts to show a tail at large 5, with clusters sizes of the 
order of a few hundreds particles. At lower temperatures, 
n{s) displays a power law regime for high values of 5, with 
a crossover point that moves to larger s with decreasing 
T. This regime is compatible with an exponent around 
—2.2 (the dashed line in the plots of Fig. 6) in agreement 
with random percolation ^ . By comparing the data at 
different volume fractions we observe that the tempera- 
ture at which the power law tail appears increases with 
increasing volume fraction. The percolation threshold, 
estimated as the temperature where 50% of the configu- 
rations percolate, monotonically increases from T ~ 0.08 
at (/) = 0.025 to T 0.14 at (/) = 0.075. This is coherent 
with the information obtained from c(n), in Fig. [4j the 
relative higher amount of bridging between chains upon 
increasing volume fraction corresponds to a higher prob- 
ability for the aggregating structure to percolate at the 
same temperature. Finally, we also conclude from Fig. 5 
that, once a percolating cluster is formed, the particles 
rapidly aggregate into a unique interconnected structure. 
This is in agreement also with the data of Fig. [4] show- 
ing that, at the lowest temperature, there are practically 
no free particles and the fraction of particles with only 
one neighbor is negligible, i.e. limited to dangling ends. 
This feature seems qualitatively the same here and in 
simulations performed on systems of smaller size (1000 
particles), suggesting that it is not much affected by the 
system size. 



D. Summary and discussion of the structural 
analysis. 

The information contained in Figs. 2-5 allows us to 
reach a convincing interpretation of the scattering pat- 
terns of Fig. 2. The behavior of c(n) indicates that the 
mesoscopic length scale discussed above corresponds to 
the formation of semi-flexible chains whose persistence 
length is directly related to the features of the effective 
interactions. The inter-molecular regime, strongly de- 
pendent on the volume fraction, corresponds instead to 
the formation of a network structure due to bridging be- 
tween different chains. Depending on the the volume 
fraction, the extended percolating structure will have dif- 
ferent amount of nodes (i.e. bridging connections) and 
therefore, presumably, significantly different mechanical 
properties. Interestingly, this corresponds to very differ- 
ent patterns in S{q) at low wave vectors: from the display 
of a wide crossover for the softest structure, to the pres- 
ence of a well defined pre-peak (i.e. a peak at wave vec- 
tors corresponding to distances significantly larger than 
the particle diameter). Differently from the assumption 



often made that such peak of correlation corresponds to 
the presence of stable finite clusters, here this pattern in- 
dicates instead the persistence length of our elongated ag- 
gregates. Our interpretation suggests that, in a case like 
this, the position of the peak could be therefore related 
to the mechanical properties of the network. Finally the 
aggregation process seems to be characterized by the fact 
that, once that a percolated structure is formed, parti- 
cles aggregate rapidly into a single network structure and 
small clusters as well as free particles completely disap- 
pear. This feature, which might be distinctive of aggre- 
gation induced colloidal gelation, significantly modifies 
the cluster size distribution and has some consequences 
also on the dynamics. 

The general understanding of the structure formation 
discussed here will be relevant to the analysis of the dy- 
namics, which is performed in the following section. 



IV. DYNAMICS 

We characterize the dynamics in our model by means 
of different quantities: The mean squared displacement 
of the particles. 



(Ar^W) = ^(E(r,(t)-r,(0))^), 



(5) 



and the incoherent scattering function 
1 ^ 

Fsiq.t) = -E(exp[iq- (r,(t) -r,(0))]). (6) 

In order to fully understand the role of the structure for- 
mation in the dynamics we also calculate time correla- 
tions of bonds and nodes, defined above, in the following 
way: 

where nij{t) = 1 if particles i and j are linked at time 
t and nij{t) = otherwise. For the nodes (i.e. particles 
connected to three other particles) Czb{t) is defined as 



Ei [{nh)-{nuY] 



(8) 



where n^i{t) = 1 if particle i is a node at time t, n^i{t) = 
otherwise. 



A. Characteristic time scales from time 
autocorrelation functions 

From the time correlation functions introduced above 
we calculate the characteristic times Ts{q) = J Fs{q^t)dt^ 
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FIG. 6: Arrhenius plot of the relaxation time Ts{q,T) 
as determined from the self-intermediate scattering function 
Fs{q,t). The solid line is a fit to the high T data for q = 3.0 
of the form r — const. /a/T- The open symbols correspond 
to different wave- vectors. The dark circles and squares are Tb 
and T36, respectively. 



Tb = J Ci}{t)dt and rsb = / Csb{t)dt^ where suitable 
cut-offs have been used to evaluate the time integrals. 
Fig. [6] shows the dependence of these different character- 
istic times on the inverse temperature at volume fraction 
(j) = 0.025 (a) and (j) = 0.075 (b). At high T the re- 
laxation dynamics at large q can be approximated by 
the function Fs{q^t) = exp[— Tg^t^/(2m)], which gives 
Tg = {^/7^n/2T)/q (solid line). At these temperatures, 
bonds are not persistent enough to create long living 
structures. Upon lowering the temperature, bond life- 
time becomes the longest relaxation time scale. More- 
over, at the lowest temperatures the lifetime of the net- 
work nodes becomes comparable to the relaxation times 
at low wave vectors. That is, the disordered network 
structure is persistent enough to affect the dynamics at 



large length scales. From the figure we also conclude that, 
whereas the bond lifetime is fairly insensitive to the vol- 
ume fraction, the node lifetime displays a stronger depen- 
dence on (j): the persistence of the nodes apparently de- 
pends on the structural features of the network, which in 
turn, as seen above, depends on (j). This already suggests, 
as discussed more in the following, that the network for- 
mation is related to the onset of cooperative dynamical 
processes. We also observe that at ^ = 0.075, as also 
found at ^ = 0.05 [25 , the relaxation time rs(g,T) at 
the smallest wave vector qmin — 0.16 shows a tempera- 
ture dependence which becomes stronger than Arrhenius 
at low T, whereas this does not happen at the lowest 
volume fraction (j) = 0.025. Upon increasing ^, the in- 
crease of node lifetime is apparently stronger. Again, this 
suggests that, whereas the bond breaking is an activated 
process mainly controlled by the potential energy param- 
eters, the node lifetime is actually affected by the features 
of the different structure at different volume fractions. 
From the information gained from S{q) and our analysis 
of section [llll we expect the large lengths scale properties 
of the structure to play a relevant role in the node life- 
time. We will come back to these considerations in the 
following, when analyzing the decay of time correlations 
of particle displacement, bonds, and nodes. 



B. Mean squared displacement 

At high temperatures, the time dependence of the par- 
ticle mean squared displacement (MSD) crosses over from 
the ballistic t^-dependence at short times to the diffusive, 
i.e. linear t-dependence. Upon lowering the temperature 
we observe the onset of a more complex behavior at a 
temperature T 0.1, at which the attraction starts to 
create persistent bonding. In Fig. [t] we plot (Ar^(t))/t 
at ^ = 0.075 for the different temperatures (main frame). 
This plot clearly shows that, upon lowering the temper- 
ature, two different localization processes arise. The first 
one, at t ^ 1, corresponds to a localization length ^ 0.2, 
which shows a relatively weak dependence on further low- 
ering temperature: this is the onset of the caging regime 
in which a particle is temporarily trapped by its near- 
est neighbors, due to the formation of bonds [25ti27] . At 
the lowest temperatures, the localization process arising 
at times t ~ 10^ corresponds to a localization distance 
of the order of 10 particles diameter, i.e. comparable 
to the size of the mesh in the disordered network struc- 
ture. In the inset, the data at the lowest temperature 
T = 0.05, where practically all the particles belong to 
a single percolated cluster, are directly compared to the 
contribution coming from the nodes of the network (par- 
ticles with coordination number n = 3), from particles 
belonging to chains (n = 2) and from the dangling ends 
(n = 1). This n-dependence gives two main indications: 
the first localization process strongly limits the mobility 
of the network nodes, but the overall MSD is dominated 
by the chain mobility; the second localization process. 
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t 



FIG. 7: Main frame: {Ar'^{t))/t as a function of time at 
volume fractions (f) — 0.075 for different temperatures. Inset: 
comparison of the data at T = 0.05 of the main frame with the 
same quantity (fuh Une) calculated for nodes, chains particles 
and dangling ends. 




t 



FIG. 8: Main frame: Particle MSD as a function of time at 
T = 0.05. Inset: MSD for the nodes at T = 0.05 for different 
volume fractions. 



present only at the lowest temperature where the life- 
time of the network nodes becomes comparable to the 
longest characteristic time scale, strongly limits the large 
scale mobility, dominating the MSD. The particle MSD 
for the different (j) at T = 0.05, i.e. when the network 
is fully developed, is plotted in Figjs] (main frame) as a 
function of time. Here it is clear that, whereas the first 
localization process (t ^ 1.) is very weakly dependent on 
0, the volume fraction dependence of the second localiza- 
tion process (t ~ 10^) is much stronger. This feature is 
also present if one isolates the contribution coming from 
the nodes of the network (inset of Figj8| at the different 
volume fractions. These observations are coherent with 
the volume fraction dependence of the structure on dif- 



ferent length scales, as shown by the S{q) in section III 



To get a deeper understanding of the correlated par- 
ticle motion, it is useful to monitor the deviation from 
a Gaussian distribution for the particle displacements, 
which is quantified by the non- Gaussian parameter 



a2{t) 



5(Ar2(t))2 



1. 



(9) 



In Fig. [9] we show 0^2 (t) as a function of time at dif- 
ferent temperatures and volume fraction ^ = 0.05: at 
high temperature there is a small peak in 0^2 (^) aris- 
ing at the crossover between the ballistic and the dif- 
fusive regime. Upon decreasing temperature, this peak 
strongly increases, indicating that the first localization 
process induces increasing non-Gaussian contributions to 
the crossover into the diffusive regime. This is actually 
similar to what is typically observed in the non- Gaussian 
parameter for supercooled liquids at the onset of the 
caging regime. At even lower temperature, at which the 
network starts to form, one can clearly see a qualitative 




FIG. 9: Non-Gaussian parameter a2{t) as a function of time 
for = 0.05 at different temperatures. 



change in the non- Gaussian contribution to particle mo- 
tion: the glassy peak intensity at t ~ 10 decreases and 
the maximum of a2{t) apparently moves towards much 
longer times, occurring approximately after the second 
localization process. This indicates that, once the aggre- 
gation leads towards the network, the main non- Gaussian 
contribution to particle motion is due to the second lo- 
calization process, i.e. on length scales of the order of 
network mesh size. This indicates that there is in fact a 
different glassy regime of the relaxation dynamics which 
is set in by the formation of the persistent network and 
whose onset is therefore marked by the second localiza- 
tion process. This observation is consistent with the find- 
ings of a recent numerical study of dipolar colloidal gels 
^40] . although in that study the role of bonds and net- 
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FIG. 10: Main frame: Non-Gaussian parameter 0^2 (t) as a 
function of time for (j) — 0.05 at T = 0.05. Inset: a2{t) at 
T = 0.05 and volume fractions = 0.025 (full line), (j) = 0.05 
(dashed line) and (j) = 0.075 (dash-dotted line). 
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FIG. 11: Fs{q,t) as a function of time at T = 1.0 and 
= 0.075 for wave vectors q = 10.0, 7.0, 5.0, 3.0, 2.0, 1.0, 
0.88, 0.54, 0.32, 0.16. The lines correspond to fitted decays 



work junctions has not been quantitatively investigated. 
We suggest that this type of glassy dynamics induced by 
the persistent network might be a distinctive feature of 
colloidal gelation. 

It is also useful to distinguish the contribution to a2 (t) 
coming from different part of the network structure as we 
do in Fig. 10 In the main frame we show a2{t)^ as cal- 
culated from all the particles at = 0.05 and T = 0.05, 
and compare it to the contribution coming from chains 
particles (n = 2) and nodes (n = 3). In spite of the fact 
that the different parts of the structure have been shown 
to give very different contributions to the dynamics, the 
data indicates that, at this low temperature, the 0^2 (t) for 
the different connectivities are quite similar. This find- 
ing is fully consistent with a recent experimental analysis 
of dynamical heterogeneities in colloidal gels of Ref. [23^ , 
where it has been proposed as an indication that dynami- 
cal heterogeneity in presence of a persistent network does 
not originate from the heterogeneous structure but is in- 
stead due to the presence of cooperative processes as in 
dense glassy systems [41 . Here we have been able to 
show that the mechanism producing this glassy regime 
is the formation of the persistent, open network. In the 
following sections we will better elucidate the nature of 
the cooperative processes. For this we use F(g, t) to in- 
vestigate the time correlations of particle motions over 
different length scales in order to understand the connec- 
tion between structural and dynamical heterogeneities in 
this system. 



C. Incoherent scattering function 



In Figs. 11 and [12] the incoherent scattering function 



Fs{q^t), as defined in Eq. (Ig]), is plotted as a function 
of the rescaled time t/rs{q)foT different wave vectors 



at = 0.075, and at T = 1.0 and T = 0.05 respec- 
tively. In agreement with the analysis performed in 
Ref. [26 at ^ = 0.05, here we also observe that at high 
temperature for q > 1.0 all the curves follow a decay 
(X exp(— (t/r(g))^) with f3 = 2 and for q < 1.0 they cross 
over to a decay with /3 = 1. This simply illustrates the 
crossover from the ballistic to the diffusive regime, tak- 
ing place on a length scale of the order of the mean free 
path and is coherent with the information obtained from 
Figs. |6j [7| and |8] At low temperature instead. Fig. 12 
shows a sharper crossover, as a function of q^ from a de- 
cay with /3 1.4 for g > 1.0 to a decay with /3 o:^ 0.55 for 
q < 1.0 (see also Ref. [^E]). At this T, finite clusters 
and free particles are rare and the main contribution to 
the displacement of the particles comes from the parti- 
cles of the network which are not directly connected to 
the nodes. This indicates that, for wave vectors q cor- 
responding to a few inter-particle distances and smaller, 
the motion of the branches of the percolating cluster is 
the main relaxation mechanism. For intermediate and 
large length scales the system shows a relaxation dynam- 
ics similar to the one found in dense glass-forming liquids, 
i.e. the presence of the disordered structure on this length 
scales leads to heterogeneous dynamics characterized by a 
stretched exponential decay of time correlations |26l |27] . 

In Fig. 13 we plot • qVT as a function of q at 
4> = 0.075 for different temperatures. The horizontal 
line in the figure, which approximates well the data for 
high T and high q^ corresponds to Tg = {^/7rm/2T)/q. 
At the lowest T, the data for large q appears to follow a 
different nea W^-ballistic regime. For T < 0.05 more than 
97% of the particles belong to one percolating cluster 
(see Fig. 5), therefore this fast regime is not really bal- 
listic and is instead due to the motion of the branches of 
the gel network [42 . The length scale q 1.0 marks the 
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FIG. 12: Fs{q,t) as a function of time at T = 0.05 and piG. 14: Fs{q,t) as a function of rescaled time at T = 1. 

= 0.075 for wave vectors q = 10.0, 7.0, 5.0, 3.0, 2.0, 1.0, (i^set) and T = 0.05 (main frame) for q = 5.0 and dif- 

0.88,0.54,0.32,0 16. The open symbols correspond to decays fg^ent volume fractions. The lines correspond to decays 

oc exp(-(t/r(^))^). ^ eM-(t/r{q)f). 



n 1 I I I I 




FIG. 13: rs{q) • qT^'^ as a function of the wave vector q at 
different temperatures and = 0.075 to illustrate the length 
scale dependence of the relaxation processes, at high and low 
T. The horizontal line corresponds to Tsq\/~{T) — ^nm/2T. 



at high (inset) and low (main frame) temperatures for 
g = 5.0. The plots indicate the change in the time decay 
from high to low T, but the ^-dependence is very weak. 
Although this is at the end not surprising in the low vol- 
ume fraction regime chosen, Fig. [15] clearly shows that 
this is not obvious. There, Fs{q^t) is plotted at high (in- 
set) and low (main frame) temperatures for q = 0.3. On 
these length scales, one can detect not only the change 
in the time decay from high to low T, but also a change 
towards a rather strong ^-dependence at low T as com- 
pared to high T. Thus, the emerging picture is basi- 
cally the dynamical analogue of Fig. [3j indicating that, 
at small length scales, i.e. up to distances of the order 
of 3-4 particle diameters, the structural heterogeneity is 
completely dominated by the interaction potential and 
therefore weakly dependent on (j): this corresponds to 
a weakly ^-dependent relaxation dynamics. At larger 
length scales instead, beyond the mesh size of the net- 
work, even relatively small changes of (j) strongly affect 
spatial correlations and the structure of the network, pro- 
ducing a strongly ^-dependent relaxation dynamics. 



crossover to a different dynamic regime. Our data shows 
that at T = 0.05 the relaxation time extracted from 
Fs{q^t) displays a strong g'— dependence, which resem- 
bles the one found in dense glass-forming liquids. How- 
ever, here this dependence is observed at wave- vectors 
that correspond not to an inter-particle distance, but to 
the mesh size of the network. 

The complex structure of the network and its strong 
heterogeneity has also a detectable effect on the (j)- 
dependence of Fs{q,t) at different wave- vectors. In 
Fig. 14, Fs{q^t) is plotted as a function of time for the 
different volume fractions (j) = 0.025, 0.05, and 0.075 



D. Bond and node time correlations 

The complex relaxation behavior just described has 
also its correspondence in the decay of time correlation 
functions for bonds and nodes of the network, as elu- 
cidated in Figs. 16 and pTj Here, Cb(t) and Csb{t)i as 
defined respectively in Eqs. ([t]) and (|8|, are plotted as 
a function of time at volume fraction <p = 0.075. The 
long time decay of Ci,{t) is well described by a single 
exponential law at all temperatures, in agreement with 
the Arrhenius dependence of the bond lifetime shown in 
Fig. [6j The activation energy related to bond breaking 
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FIG. 15: Fs{q,t) as a function of time at T = 1.0 (inset) 
and T — 0.05 (main frame) for q — 0.3 and different volume 
fractions. 




FIG. 16: The bond time correlation function Chit) as a func- 
tion of time, at volume fraction <j) — 0.075. T — 1.0, 0.5, 0.2, 
0.1, 0.09, 0.08, 0.07, 0.06, and 0.05 from left to right. The long 
time decay follows an exponential law at all temperatures. 



does not significantly depend on the volume fraction, as 
expected. As already noticed in Ref. |27], we can actually 
distinguish two different regimes, one at high tempera- 
tures (T > 0.1) and a second one at low temperature 
(T < 0.1). At high temperatures the particle collisions 
promote uncorrelated bond breaking or formation, lead- 
ing to a short time decay of bond correlation and hence 
to a quick decay of correlations. At low temperatures 
instead, the role of particle collisions decreases and the 
energy activation process becomes the only relevant pro- 
cess in the decay of bond correlation. In Fig. [T7| the 
time correlation function Cs}){t) for particles with coor- 
dination number n = 3 is plotted as a function of time. 
At high temperatures (T > 0.1) particles of connectiv- 
ity 3 are rare (see Fig. H|, and therefore the statistics is 




FIG. 17: The time correlation function Cstit) as a function 
of time, at volume fraction = 0.075. T = 1.0, 0.5, 0.2, 
0.1, 0.09, 0.08, 0.07, 0.06, and 0.05 from left to right. At the 
lowest temperatures, where it is the time correlation function 
of the network nodes, C^hif) displays a strongly stretched 
decay (/3 = 0.54). 



rather poor. At temperatures T > 0.07, the long time 
decay of time correlation functions Csb{t) follows a sim- 
ple exponential law, with a characteristic relaxation time 
rsi){(j),T) increasing with decreasing T. This is coherent 
with the Arrhenius dependence of r^b in Figj6] In con- 
trast to this, at temperatures at which a persistent span- 
ning network is present in the system, i.e. T < 0.06, the 
decay of Csbit) becomes stretched, with a stretching ex- 
ponent P which decreases with T {P 0.54 at T = 0.05). 
This indicates that, once the network is formed and it is 
sufficiently persistent, the process of breaking and forma- 
tion of the nodes is associated not only to the overcoming 
of an activation energy but also to a heterogeneous and 
cooperative dynamic process [25] . 

We have further investigated the type of cooperative 
processes in which the nodes might be involved, in par- 
ticular by measuring the occurrence of events at which 
the breaking of one node makes that one of its neigh- 
bors of connectivity 2 becomes itself a node. This would 
make the nodes to slide along the chains. We have found 
that in the network regime at least 25% of the events 
corresponds to this situation. This behavior points out 
a further similarity, in spite of the very different type of 
interactions, between the system studied here and the 
dipolar colloidal gels of Ref. [40 , in which one could eas- 
ily expect these events to occur. These results elucidate 
well the complex interplay between structure and dynam- 
ics which is fundamental in these systems: Although the 
slow dynamics is due to the bond lifetime becoming suf- 
ficiently long (see Figj6|, it is only the combination of 
sufficiently persistent bonds together with the branching 
of the aggregates, i.e. the formation of a stress bear- 
ing network, which eventually determines the arising of 
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large scale cooperative relaxation processes at these low 
volume fractions. 



E. Summary and discussion of the dynamical 
properties 

We can now put together all the elements obtained 
from the different quantities into a unique coherent 
picture. The onset of the aggregation, signaled for 
example by the change of shape of the cluster size 
distribution around T = 0.1 in Fig. 5, corresponds to 
the onset of a dynamical regime at which relaxation 
processes are dominated by bond-breaking as opposed 
to particle collisions (Fig. [6|. This dynamical regime 
may resemble, to some extent, the onset of caging in 
dense systems: The MSD, signals a localization process 
over similar length scales (Figs, [t] and |8| accompanied 
by a growing degree of heterogeneity in the distribution 
of particle displacements (Figs, [o] and [lO|, as compared 
to a Gaussian one. In spite of these similarities, we have 
shown that here this regime does not imply structural 
arrest. The heterogeneity of the dynamics is clearly 
arising from the wide distribution of the sizes of the 
aggregates [^61 [43ti45] . Upon further decreasing tem- 
perature, the aggregation process leads to the formation 
of an interconnected network. Once that the network 
is persistent enough, i.e. the node lifetime starts to 
be comparable to the longest relaxation times in the 
system (Fig. [6|, a new, slow dynamic regime sets in. 
This has certainly the hallmark of gelation, due to the 
formation of the locally rigid, persistent network, and it 
is also indicated by the strong localization process in the 
MSD, strikingly dominated by the network nodes (inset 
of Fig. [7|), over length scales of the order of the mesh 
size of the network. This second dynamical regime is 
associated to a qualitative change in dynamical hetero- 
geneity (Fig. [9|. The contribution arising from the first 
localization process (caging) becomes negligible: there 
is a new significant non- Gaussian contribution to the 
distribution of particle displacements arising from the 
second localization process (see the curve at T = 0.07 
in Fig. |9| and moving towards longer and longer time 
scales. The heterogeneity of this slow dynamics has no 
straightforward relation to the structure (Fig.[lO|, there- 
fore suggesting that it signals instead the presence of 
new cooperative processes. The study of the relaxation 
dynamics over different length scales elucidates well that 
this second dynamical regime is characterized by the 
coexistence of very different relaxation processes over 
different length scales (Figs. TT][T5 ): fast motion of pieces 
(chains and dangling ends) of the gel at small distances 
and slow, stretched exponential processes related to the 
network rearrangements at length scales larger than the 



network mesh size. This clearly indicates the network 
origin of this dynamical regime, which, being character- 
ized by slower and slower, complex relaxation, points 
to structural arrest. Our comparative analysis indicates 
therefore that this complex dynamics has glassy features 
and strong similarities to the one of dense systems, with 
the difference that the strong coupling in particle motion 
(which originates the structural arrest) is not induced 
by the crowding but by the presence of the persistent 
network. The fact that the onset of very slow, stretched 
exponential relaxations only takes place at sufficiently 
low wave vectors is a consequence of that. Finally, by 
looking to bond and nodes relaxation we have found 
that the cooperative processes underlying this network 
induced glassy dynamics is intimately related to the 
network nodes (Fig. 



V. CONCLUSIONS 

We have discussed the behavior of a model for col- 
loidal gels based on the presence of directional effective 
interactions. In this model, the aggregation leads to an 
open persistent network structure without imposing a 
fixed connectivity to the gel units. With these premises, 
our study gives new insights into the physics of colloidal 
gelation. From the structural point of view, we have 
shown that the aggregation process takes place via a 
random percolation mechanism, but once a percolating 
structure is formed, it rapidly evolves towards a persis- 
tent, fully connected open network. Also in this case, 
gelation seems associated to the presence of a pre-peak 
in the static structure factor, here directly related to the 
network structure. Our comparative analysis of struc- 
ture and dynamics indicates that the persistent network 
introduces slow, cooperative processes intimately related 
to the network nodes and sets in a peculiar kind of glassy 
dynamics. This scenario is also consistent with the re- 
sults of a recent numerical study of diluted dipolar col- 
loidal gels . We think that this points to the formation 
of a persistent network as the mechanism responsible for 
the onset of the glassy dynamics in colloidal gels and ex- 
plains the close connection between gelation and glassy 
structural arrest typically observed in these systems. We 
therefore suggest that this scenario may in fact be rel- 
evant to the physics of colloidal gels on a more general 
ground. 
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